Description
Uses a collapsed (marginalizes over a, b,
amp numerically) gibbs sampler which samples from the
conditionals with slice sampling. Also re-parameterizes density with
(sqrt(gwidth^2+lwidth^2), gwidth-lwidth) instead of
(gwidth, lwitdh) if
min(gwidth/lwidth, lwidth/gwidth) > 0.05 (otherwise
transformation becomes unstable). Finds global maxima in posterior using
population based golden section search and uses it as the starting point
for the slice sampler. About x10 times faster than the previous
implementation. Using Zell Prior could be another x10.
Data
library(spectralmc)
#toluene data
cdata <- get_toluene_measured_data()
x <- cdata$x
y <- cdata$y
Algorithm
peak_spawner <- function() {
params <- list(
amp = rnorm(1, 1, 1), # doesnt matter
lwidth = rnorm(1, 0.1, 0.001), # matters!!!
gwidth = rnorm(1, 5, 0.001), # matters!!!
pos = runif(1, 1, 1), # doesnt matter
a = 0, # doesnt matter
b = 0 # doesnt matter
)
return(params)
}
fit_kp <- 30
max_a <- 1e-2
lower <- c(-Inf, -max_a, 0)
upper <- c(Inf, max_a, Inf)
noise_sigma <- 100 # estimated by eye; intentionally to big such that tiny peaks are overnoised
signal_model <- voigt.model
res <- chains.rjmcmc_like_collapsed_slicer_gibbs(
x,
y,
signal_model,
peak_spawner,
lower,
upper,
noise_sigma,
fit_kp = fit_kp,
add_peak_every = 10,
iter = 20000,
print_bar = TRUE,
half_steps = TRUE,
decorr_trafo = TRUE,
save_loc = save_loc,
checkpoint_every = 50,
)
Results
Fit
All
plt.fit.chain.interactive(x, y, res$samples, signal_model, step_width = 200)
Burn-in
plt.fit.chain.interactive(x, y, res$samples[seq(1, 500)], signal_model, step_width = 5)
Traceplots
plt.traceplot(res$samples, "pos")
plt.traceplot(res$samples, "amp")
plt.traceplot(res$samples, "lwidth")
plt.traceplot(res$samples, "gwidth")
plt.traceplot(res$samples, "a")
plt.traceplot(res$samples, "b")
Posterior plot
plt.posterior(res$post_vals, steps = NULL)
Todo
- Critical Bug amp, a, b are not sampled instead they
are set to roughly their means (not really either); Easy to fix.
- Not a bug: Uses a safeguard which prevents messing
up the fit if the slice sampler fails, this however destroys the gibbs
sampler property and causes rejection plateaus.
- Done: nu, be transformation becomes bad (why?) if
either lwidth or gwidth is really small compared to the other, should
switch to nu,be only if
min(lwidth/gwidth, gwidth/lwidth) > 0.05
LS0tDQp0aXRsZTogIkNvbGxhcHNlZCBHaWJicyBTYW1wbGVyIHdpdGggR1NTIFNsaWNlIFNhbXBsaW5nIg0KYXV0aG9yOiBKYW4gTWVpw59uZXJeW1JXVEggQWFjaGVuLCBwaGlsaXBwLm1laXNzbmVyQHJ3dGgtYWFjaGVuLmRlXQ0Kb3V0cHV0Og0KICBodG1sX25vdGVib29rOg0KICAgIHRvYzogdHJ1ZQ0KICAgIHRvY19kZXB0aDogNQ0KICAgIHRvY19mbG9hdDoNCiAgICAgIGNvbGxhcHNlZDogZmFsc2UNCiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUNCi0tLQ0KDQo8c2NyaXB0IHR5cGU9InRleHQvamF2YXNjcmlwdCI+DQogIGZ1bmN0aW9uIGp1bXBfaGVhZGVyKGtleSl7DQogICAgJCgnOmhlYWRlcjpjb250YWlucygnK2tleSsnKScpWzBdLnNjcm9sbEludG9WaWV3KCk7DQogIH0NCg0KICBmdW5jdGlvbiBqdW1wX21hcmtlZChrZXkpew0KICAgICQoJyMnICsga2V5KVswXS5zY3JvbGxJbnRvVmlldygpOw0KICB9DQoNCiAgZnVuY3Rpb24gcGxvdFpvb20oZWwpew0KICAgICAgaWYoZG9jdW1lbnQuZnVsbHNjcmVlbikgew0KICAgICAgICBkb2N1bWVudC5leGl0RnVsbHNjcmVlbigpDQogICAgICB9IGVsc2Ugew0KICAgICAgICAkKGVsKS5jbG9zZXN0KCcuanMtcGxvdGx5LXBsb3QnKVswXS5yZXF1ZXN0RnVsbHNjcmVlbigpOw0KICAgICAgfQ0KICB9DQoNCiAgJCggZG9jdW1lbnQgKS5yZWFkeShmdW5jdGlvbigpIHsNCiAgICAkKCIubW9kZWJhci1idG4ucGxvdGx5anNpY29uLm1vZGViYXItYnRuLS1sb2dvIikucmVwbGFjZVdpdGgoDQogICAgYA0KICAgIDxhIHJlbD0idG9vbHRpcCIgb25jbGljaz1wbG90Wm9vbSh0aGlzKSBjbGFzcz0ibW9kZWJhci1idG4gZnVsbHNjcmVlbi1idG4iIGRhdGEtdGl0bGU9IkZ1bGwgU2NyZWVuIiBkYXRhLWF0dHI9Inpvb20iIGRhdGEtdmFsPSJhdXRvIiBkYXRhLXRvZ2dsZT0iZmFsc2UiIGRhdGEtZ3Jhdml0eT0ibiIgPg0KICAgICAgPHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHZpZXdCb3g9IjAgMCA0NDggNTEyIiBjbGFzcz0iaWNvbiIgaGVpZ2h0PSIxZW0iIHdpZHRoPSIxZW0iPg0KICAgICAgICA8IS0tISBGb250IEF3ZXNvbWUgUHJvIDYuMS4xIGJ5IEBmb250YXdlc29tZSAtIGh0dHBzOi8vZm9udGF3ZXNvbWUuY29tIExpY2Vuc2UgLSBodHRwczovL2ZvbnRhd2Vzb21lLmNvbS9saWNlbnNlIChDb21tZXJjaWFsIExpY2Vuc2UpIENvcHlyaWdodCAyMDIyIEZvbnRpY29ucywgSW5jLiAtLT4NCiAgICAgICAgPHBhdGggZD0iTTEyOCAzMkgzMkMxNC4zMSAzMiAwIDQ2LjMxIDAgNjR2OTZjMCAxNy42OSAxNC4zMSAzMiAzMiAzMnMzMi0xNC4zMSAzMi0zMlY5Nmg2NGMxNy42OSAwIDMyLTE0LjMxIDMyLTMyUzE0NS43IDMyIDEyOCAzMnpNNDE2IDMyaC05NmMtMTcuNjkgMC0zMiAxNC4zMS0zMiAzMnMxNC4zMSAzMiAzMiAzMmg2NHY2NGMwIDE3LjY5IDE0LjMxIDMyIDMyIDMyczMyLTE0LjMxIDMyLTMyVjY0QzQ0OCA0Ni4zMSA0MzMuNyAzMiA0MTYgMzJ6TTEyOCA0MTZINjR2LTY0YzAtMTcuNjktMTQuMzEtMzItMzItMzJzLTMyIDE0LjMxLTMyIDMydjk2YzAgMTcuNjkgMTQuMzEgMzIgMzIgMzJoOTZjMTcuNjkgMCAzMi0xNC4zMSAzMi0zMlMxNDUuNyA0MTYgMTI4IDQxNnpNNDE2IDMyMGMtMTcuNjkgMC0zMiAxNC4zMS0zMiAzMnY2NGgtNjRjLTE3LjY5IDAtMzIgMTQuMzEtMzIgMzJzMTQuMzEgMzIgMzIgMzJoOTZjMTcuNjkgMCAzMi0xNC4zMSAzMi0zMnYtOTZDNDQ4IDMzNC4zIDQzMy43IDMyMCA0MTYgMzIweiIvPg0KICAgICAgPC9zdmc+DQogICAgPC9hPg0KICAgIGApOw0KICB9KTsNCjwvc2NyaXB0Pg0KDQoNCiMgRGVzY3JpcHRpb24NClVzZXMgYSBjb2xsYXBzZWQgKG1hcmdpbmFsaXplcyBvdmVyIGBhYCwgYGJgLCBgYW1wYCBudW1lcmljYWxseSkNCmdpYmJzIHNhbXBsZXIgd2hpY2ggc2FtcGxlcyBmcm9tIHRoZSBjb25kaXRpb25hbHMgd2l0aCBzbGljZSBzYW1wbGluZy4gQWxzbyByZS1wYXJhbWV0ZXJpemVzIGRlbnNpdHkgd2l0aCBgKHNxcnQoZ3dpZHRoXjIrbHdpZHRoXjIpLCBnd2lkdGgtbHdpZHRoKWAgaW5zdGVhZCBvZiBgKGd3aWR0aCwgbHdpdGRoKWAgaWYgYG1pbihnd2lkdGgvbHdpZHRoLCBsd2lkdGgvZ3dpZHRoKSA+IDAuMDVgIChvdGhlcndpc2UgdHJhbnNmb3JtYXRpb24gYmVjb21lcyB1bnN0YWJsZSkuIEZpbmRzIGdsb2JhbCBtYXhpbWEgaW4gcG9zdGVyaW9yIHVzaW5nIHBvcHVsYXRpb24gYmFzZWQgZ29sZGVuIHNlY3Rpb24gc2VhcmNoIGFuZCB1c2VzIGl0IGFzIHRoZSBzdGFydGluZyBwb2ludCBmb3IgdGhlIHNsaWNlIHNhbXBsZXIuIEFib3V0IHgxMCB0aW1lcyBmYXN0ZXIgdGhhbiB0aGUgcHJldmlvdXMgaW1wbGVtZW50YXRpb24uIFVzaW5nIFplbGwgUHJpb3IgY291bGQgYmUgYW5vdGhlciB4MTAuDQoNCiMgRGF0YQ0KYGBge3IgZWNobyA9IEZBTFNFfQ0KI3NldHdkKCdDOi9Vc2Vycy9KYW4vRGVza3RvcC9JU1dfU3BlY3RyYUJheWVzL0lTV19TcGVjdHJhQmF5ZXMvYmF5ZXNfbWNtYycpDQoNCmdldF90b2x1ZW5lX21lYXN1cmVkX2RhdGEgPC0gZnVuY3Rpb24oKXsNCiAgcGF0aCA8LSAiQzovVXNlcnMvSmFuL0Rlc2t0b3AvSVNXX1NwZWN0cmFCYXllcy9JU1dfU3BlY3RyYUJheWVzL2JheWVzX21jbWMvbm90ZWJvb2tzL2dvbGRlbl9naWJic190b2x1bmUvdG9sdWVuZV9tZWFzdXJlZC5jc3YiDQogIHRvbHVlbmVfbWVhc3VyZWQgPC0gcmVhZC5jc3YocGF0aCwgaGVhZGVyPUZBTFNFLCBzZXA9IjsiKQ0KICB4IDwtIGFzLnZlY3Rvcih0b2x1ZW5lX21lYXN1cmVkWywxXSkNCiAgeXRydWUgPC0gYXMudmVjdG9yKHRvbHVlbmVfbWVhc3VyZWRbLDJdKQ0KICByZXR1cm4obGlzdCh4PXgseT15dHJ1ZSkpDQp9DQpgYGANCmBgYHtyfQ0KbGlicmFyeShzcGVjdHJhbG1jKQ0KDQojdG9sdWVuZSBkYXRhDQpjZGF0YSA8LSBnZXRfdG9sdWVuZV9tZWFzdXJlZF9kYXRhKCkNCnggPC0gY2RhdGEkeA0KeSA8LSBjZGF0YSR5DQpgYGANCiMgQWxnb3JpdGhtDQpgYGB7ciBldmFsID0gRkFMU0V9DQpwZWFrX3NwYXduZXIgPC0gZnVuY3Rpb24oKSB7DQogIHBhcmFtcyA8LSBsaXN0KA0KICAgICAgICAgIGFtcCA9IHJub3JtKDEsIDEsIDEpLCAjIGRvZXNudCBtYXR0ZXINCiAgICAgICAgICBsd2lkdGggPSBybm9ybSgxLCAwLjEsIDAuMDAxKSwgIyBtYXR0ZXJzISEhDQogICAgICAgICAgZ3dpZHRoID0gcm5vcm0oMSwgNSwgMC4wMDEpLCAjIG1hdHRlcnMhISENCiAgICAgICAgICBwb3MgPSBydW5pZigxLCAxLCAxKSwgIyBkb2VzbnQgbWF0dGVyDQogICAgICAgICAgYSA9IDAsICMgZG9lc250IG1hdHRlcg0KICAgICAgICAgIGIgPSAwICMgZG9lc250IG1hdHRlcg0KICApDQogIHJldHVybihwYXJhbXMpDQp9DQoNCmZpdF9rcCA8LSAzMA0KbWF4X2EgPC0gMWUtMg0KbG93ZXIgPC0gYygtSW5mLCAtbWF4X2EsIDApDQp1cHBlciA8LSBjKEluZiwgbWF4X2EsIEluZikNCm5vaXNlX3NpZ21hIDwtIDEwMCAjIGVzdGltYXRlZCBieSBleWU7IGludGVudGlvbmFsbHkgdG8gYmlnIHN1Y2ggdGhhdCB0aW55IHBlYWtzIGFyZSBvdmVybm9pc2VkDQpzaWduYWxfbW9kZWwgPC0gdm9pZ3QubW9kZWwNCg0KcmVzIDwtIGNoYWlucy5yam1jbWNfbGlrZV9jb2xsYXBzZWRfc2xpY2VyX2dpYmJzKA0KICAgICAgICB4LA0KICAgICAgICB5LA0KICAgICAgICBzaWduYWxfbW9kZWwsDQogICAgICAgIHBlYWtfc3Bhd25lciwNCiAgICAgICAgbG93ZXIsDQogICAgICAgIHVwcGVyLA0KICAgICAgICBub2lzZV9zaWdtYSwNCiAgICAgICAgZml0X2twID0gZml0X2twLA0KICAgICAgICBhZGRfcGVha19ldmVyeSA9IDEwLA0KICAgICAgICBpdGVyID0gMjAwMDAsDQogICAgICAgIHByaW50X2JhciA9IFRSVUUsDQogICAgICAgIGhhbGZfc3RlcHMgPSBUUlVFLA0KICAgICAgICBkZWNvcnJfdHJhZm8gPSBUUlVFLA0KICAgICAgICBzYXZlX2xvYyA9IHNhdmVfbG9jLA0KICAgICAgICBjaGVja3BvaW50X2V2ZXJ5ID0gNTAsDQopDQpgYGANCmBgYHtyICBlY2hvPUZBTFNFfQ0Kc2V0d2QoJ0M6L1VzZXJzL0phbi9EZXNrdG9wL0lTV19TcGVjdHJhQmF5ZXMvSVNXX1NwZWN0cmFCYXllcy9iYXllc19tY21jJykNCnJlcyA8LSByZWFkUkRTKCJub3RlYm9va3MvZ29sZGVuX2dpYmJzX3RvbHVuZS9yZXMuUkRhdGEiKQ0KZml0X2twIDwtIDMwDQptYXhfYSA8LSAxZS0yDQpsb3dlciA8LSBjKC1JbmYsIC1tYXhfYSwgMCkNCnVwcGVyIDwtIGMoSW5mLCBtYXhfYSwgSW5mKQ0Kbm9pc2Vfc2lnbWEgPC0gMTAwICMgZXN0aW1hdGVkIGJ5IGV5ZTsgaW50ZW50aW9uYWxseSB0byBiaWcgc3VjaCB0aGF0IHRpbnkgcGVha3MgYXJlIG92ZXJub2lzZWQNCnNpZ25hbF9tb2RlbCA8LSB2b2lndC5tb2RlbA0KYGBgDQojIFJlc3VsdHMNCg0KIyMgRml0DQoNCiMjIyBBbGwNCmBgYHtyfQ0KcGx0LmZpdC5jaGFpbi5pbnRlcmFjdGl2ZSh4LCB5LCByZXMkc2FtcGxlcywgc2lnbmFsX21vZGVsLCBzdGVwX3dpZHRoID0gMjAwKQ0KYGBgDQoNCiMjIyBCdXJuLWluDQpgYGB7cn0NCnBsdC5maXQuY2hhaW4uaW50ZXJhY3RpdmUoeCwgeSwgcmVzJHNhbXBsZXNbc2VxKDEsIDUwMCldLCBzaWduYWxfbW9kZWwsIHN0ZXBfd2lkdGggPSA1KQ0KYGBgDQojIyMgVHJhY2VwbG90cw0KYGBge3J9DQpwbHQudHJhY2VwbG90KHJlcyRzYW1wbGVzLCAicG9zIikNCnBsdC50cmFjZXBsb3QocmVzJHNhbXBsZXMsICJhbXAiKQ0KcGx0LnRyYWNlcGxvdChyZXMkc2FtcGxlcywgImx3aWR0aCIpDQpwbHQudHJhY2VwbG90KHJlcyRzYW1wbGVzLCAiZ3dpZHRoIikNCnBsdC50cmFjZXBsb3QocmVzJHNhbXBsZXMsICJhIikNCnBsdC50cmFjZXBsb3QocmVzJHNhbXBsZXMsICJiIikNCmBgYA0KIyMjIFBvc3RlcmlvciBwbG90DQpgYGB7cn0NCnBsdC5wb3N0ZXJpb3IocmVzJHBvc3RfdmFscywgc3RlcHMgPSBOVUxMKQ0KYGBgDQoNCiMgVG9kbw0KLSAqKkNyaXRpY2FsIEJ1ZyoqIGFtcCwgYSwgYiBhcmUgbm90IHNhbXBsZWQgaW5zdGVhZCB0aGV5IGFyZSBzZXQgdG8gcm91Z2hseSB0aGVpciBtZWFucyAobm90IHJlYWxseSBlaXRoZXIpOyBFYXN5IHRvIGZpeC4NCi0gKipOb3QgYSBidWcqKjogVXNlcyBhIHNhZmVndWFyZCB3aGljaCBwcmV2ZW50cyBtZXNzaW5nIHVwIHRoZSBmaXQgaWYgdGhlIHNsaWNlIHNhbXBsZXIgZmFpbHMsIHRoaXMgaG93ZXZlciBkZXN0cm95cyB0aGUgZ2liYnMgc2FtcGxlciBwcm9wZXJ0eSBhbmQgY2F1c2VzIHJlamVjdGlvbiBwbGF0ZWF1cy4NCi0gKipEb25lKio6IG51LCBiZSB0cmFuc2Zvcm1hdGlvbiBiZWNvbWVzIGJhZCAod2h5PykgaWYgZWl0aGVyIGx3aWR0aCBvciBnd2lkdGggaXMgcmVhbGx5IHNtYWxsIGNvbXBhcmVkIHRvIHRoZSBvdGhlciwgc2hvdWxkIHN3aXRjaCB0byBudSxiZSBvbmx5IGlmIGBtaW4obHdpZHRoL2d3aWR0aCwgZ3dpZHRoL2x3aWR0aCkgPiAwLjA1YA0K